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We propose a methodology to design optimal pulses for achieving quantum optimal control on molecular 
systems. Our approach constrains pulse shapes to linear combinations of a fixed number of experimentally 
relevant pulse functions. Quantum optimal control is obtained by maximizing a multi-target fitness function 
with genetic algorithms. As a first application of the methodology we generated an optimal pulse that 
successfully maximized the yield on a selected dissociation channel of a diatomic molecule. Our pulse is 
obtained as a linear combination of linearly chirped pulse functions. Data recorded along the evolution of the 
genetic algorithm contained important information regarding the interplay between radiative and diabatic 
processes. We performed a principal component analysis on these data to retrieve the most relevant processes 
along the optimal path. Our proposed methodology could be useful for performing quantum optimal control 
on more complex systems by employing a wider variety of pulse shape functions. 
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I. INTRODUCTION 

Quantum optimal control (QOC) methodologies are em¬ 
ployed to produce optimal pulses capable of steering the 
quantum state of a molecule towards a given target state. 
QOC is accomplished by using the matter-radiation 
interaction and the time-frequency coherence of lightP 2 I 

Among all theoretical QOC approaches, gradient-based 
feedback (GBF) methods are widely employed due 
to their robustness and rapid convergence? Although 
GBF approaches have been used to compute pulses 
capable of attaining desired target molecular states,®® 
their applicability still faces serious limitations. On 
one hand, the convergence of GBF approaches towards 
the closest optimum is biased by the choice of the 
initial guess for the field. On the other hand, the 
precise reproduction in the laboratory of the generated 
optimal pulses is still challenging due to their complexity. 

In the past two decades, alternative QOC approaches 
based on heuristic optimization strategies have been 
proposed to improve the convergence towards the global 
optimum. Among these heuristic schemes, Genetic Al¬ 
gorithms (GA) have become the methods of choice?®® 
For instance, optimal pulses in the frequency domain 
have been generated with the aid of GA to control 
molecular eve nts relevant for quantum information 
processing. 10 In spite of overcoming the limitation 
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of GBF methods to reach the global optimum, cur¬ 
rent GA implementations still fail at generating optimal 
pulses that can be precisely reproduced in the laboratory. 

A desirable QOC approach designed to overcome the 
aforementioned limitations must guarantee the conver¬ 
gence towards the global optimum and constrain the 
pulse generation to linear combinations of experimen¬ 
tally attainable pulse functions. 

We propose a novel QOC method designed to overcome 
these limitations. Our approach introduces a GA 
methodology to optimize the parameters of a linear 
combination of analytical pulse shapes, which are ex¬ 
perimentally synthetizable with linear optics techniques. 
As a first test of the efficiency of the method, we control 
the dissociation yield of a diatomic molecule through a 
selected channel. In addition to controlling molecular 
processes, the proposed methodology generates data 
that can be used to analyze the interplay between 
diabatic and photodynamic processes along the optimal 
path. To this aim we perform a statistical principal 
component analysis (PCA) on the distribution of time 
integrated transition moment integrals calculated during 
the GA evolution. The resulting PCA singular values 
and vectors provide key information about the intricate 
sequence of events driven by the optimal pulse. 


The remainder of this paper is organized as follows: sec¬ 
tion [0] introduces the notation, the theoretical back¬ 
ground and the proposed GA methodology; section [Til 
shows the application of this methodology to perform 
optimal control of the dissociation yields of a model di¬ 
atomic molecule; section [IV] provides some concluding 
remarks and perspectives of this work. 
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II. THEORY 


The wavefunction, T, of a molecular sys¬ 
tem composed of N nuclei with coordinates 
R = (Ri, R 2 ,..., Ri,..., Rn), and n electrons 
with coordinates r = ( 1 * 1 , r 2 ,..., n,..., r n ) can be 

calculated by solving the time-dependent Schrodinger 
equation (atomic units will be used hereafter) 


<9T 

dt 


-iH T, 


(i) 


here H = Ho + H int , where Ho is the molecular Hamil¬ 
tonian, and H int is the interaction of the system with 
the pulse. 


Under the Born-Oppenheimer ansatz^for the wavefunc¬ 
tion of Eq. [l] the evolution of the nuclear wavepacket, 
Xk , can be obtained by tracing out the electronic degrees 
of freedom: 

[T(R) + -E’fe(R)] Xfc(R; 0+E c kiXk( R; t) = i^Xk( R; t), 

(2a) 


The nuclear state of the molecule, |X(R;£)), can be ex¬ 
pressed in a diabatic basis as a vector of nuclear config¬ 
urations: 


\X(R-t)) 


( xi(R; t) 

\X-/v(R;£) 


( 5 ) 


The formal solution of Eq. [T| is the Green’s function of 
the system of coupled equations of Eq. |2j 

|X(R; t + 2 St)) = e ~ m2St |X(R; t )). (6) 


This Green’s function can be approximated- 3 in terms of 
a symmetric double Trotter expansion to obtain: 

^—iTL28t ^ g— iTStg—iV25tiTSt /y\ 


The propagation of the initial (diabatic) nuclear 
wavepacket for each St, is achieved by switching back and 
forth between the diabatic and adiabatic representations 
of the wavepacket 


A. Genetic Algorithm 


Cu( R) = <* fc | T(R) I**) ~EWj (<H| Vj l$i> Vj ’ 

(2b) 

here Efc(R) is the potential energy surface (PES) 
corresponding to electronic state |<£fc), and Cki{ R) are 
diabatic couplings between electronic states | <£>/-) and 

!*«>• 


The radiation-matter interaction, H int , is considered 
semiclassically within the dipolar approximation, 

?C‘( R ,t)=*~4:Trn k i(R)e(t), (3) 

s(t) = Y^Vi(t), 

i 


The GA scores the individuals using a fitness function, 
Jq. For an observable O the fitness function is 


J 6 = #[$] + J (2) [e(f)], 


J o 

jW[^] = ('I'(T) 6 4-(T) 


o 


jM[e(t)} = - f \\e{t)\\ 2 dt 

Jo 


(8) 

( 9 ) 

( 10 ) 


here J i2> [e(t )] penalizes the fluence of the laser. We pro¬ 
pose the multi-target fitness function as: 


k 

J = J o 1 + E j— 

i =2 J Oi 


( 11 ) 


where (iki is an element of the electronic transition 
dipole matrix. 


to maximize observable 0\, while minimizing observ¬ 
ables Oi^L\. 


We propose to express the pulse, e(t ), as a linear combi¬ 
nation of linearly chirped pulse (LCP) functions. These 
LCPs have been used successfully to control photochemi¬ 
cal reactions. 15 22 The time profile of a LCP can be writ¬ 
ten as: 


V(t) = Eq exp 


( t~rof \ 
2 t 2 ) 


COS ( -c(t - To) +Uo(t- To) 

( 4 ) 


here To is the time shift, Eo is the pulse amplitude, c is 
the chirp rate, ujo is the central frequency, and r is the 
pulse width. As observed, the instantaneous frequency 
of a LCP changes linearly with time. 


In our GA implementation, the chromosome of the 
i-th individual is built by concatenating the vectors 
of genes of the N LCPs, = { 2 / 1 , 2/27 • “ ■ , 2 /jv}- Each 
LCP is represented by a 5-vector yj = (Eo,to,c,t,ujo)- 
We propose two genetic operations to evolve the N 
LCPs: (1) a mutation to change one or more genes with 
,probability 7 and (2) a crossover operation acting 
on two individuals to generate a new individual with 
probability 7Vx • The new chromosome is generated as a 
fit ness-weighted linear combination of the parents. 

Optimizations were carried out following the next steps: 

0. Initialize the population with a list of individuals 
whose chromosomes were generated randomly. 25 
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1. Propagate individuals using the same initial con¬ 
dition and calculate their scores with the fitness 
function of Eq. EH 

2 . Sort individuals in descending order according to 
their fitness. 

3. Compute the cumulative fitness of the entire pop¬ 
ulation. 

4. Keep the individuals with fitness greater or equal 
than the cumulative fitness. Discard the remainder. 

5. Replace those individuals discarded in step 4 with 
new ones generated by crossing over the survivors 
with probability tvx- 

6 . Go through the new population list sampling the 
probability of mutation of each gene with a Monte 
Carlo scheme. Mutate those genes with probability 
lower than i tm- 

7. Go back to step 1 if the maximum number of gen¬ 
erations has not been reached. 


III. NUMERICAL TEST 

As a first application of the proposed methodology we 
selectively controlled the yields through the dissociation 
channels of a model diatomic molecule.^ This model 
could be related to the control over the yield of homolytic 
or heterolytic bond cleavageP 

The molecule is modeled with three PESs: Ei(R), E 2 (R), 
and Es(R) which are coupled diabatically™ The off- 
diagonal terms of the nuclear Hamiltonian matrix are 

7*12 =Mi 2 (R)e(t), ( 12 ) 

«13 =M !3 (R)e(t), (13) 

#23 = Kf + C 23 (R) = M23 (R)e(t) + C 23 (R). (14) 

Table [T] lists all the terms and parameters employed in 
the time propagation of the molecular wavepacket. 

We initialize the propagation by placing a Gaussian 
wavepacket on the PES Es(R), with center at R = r e 
and width 0.265. Th e free evolution of this wavepacket 
is displayed in Figure [T ] 27 Panel A shows the population 
on the i-th PES, (xi \ Xi)i while panel B displays the 
currents J 2 and J 3 , calculated using 

Ji = Jlj 0 ^ R d) dt - ( 15 ) 

Here p is the momentum operator, Rd is the dissociation 
limit distance, and T is the final time of the propagation. 
As observed in Figure [l] the norm is not conserved along 
the propagation since we imposed absorbing boundary 
conditions by employing optical potentials. An inter¬ 
esting feature of the free evolution of the population 


TABLE I. Parameters defining the PESs and diabatic cou¬ 
plings. Time and space discretization values are provided. 


Ground state potential 

E 1 (R) = D e (l - exp (—y(R — r e ))) 2 

D e = 0.075 

o; e = 0.0077782 

r e = 4.125 

M = 1836.0 

y LO e \/M/2/D e _ 

Dissociative potentials 

E ot (R) — m cx (R — 3.15) + b a + cexp (—d(R — 3.15)) 
m 2 = -0.008681 m 3 = -0.01736 
b 2 = 0.1805 b 3 = 0.2 
c = 0.25 

d = 5.0_ 

Dipole functions and diabatic couplings 
M 12 (R) — M 23 (R) = /xexp (-a(R - r x )‘ 2 ) 

/ii 3 (R) = /xtanh (—10.0(8.0 — R)) r > r x 
/jli 3 (R) = /xtanh (—10.0(r — 3.0)) r <r x 

li = 0 . 2 ; 

a = 5.0 
r x — 5.25 

C 23 (R) = A exp (~a(R - r x ) 2 ) 

Intermediate coupling: A = 0.003 
Strong coupling: A = 0.0075 

Optical potential^^ 

Grid spacing Ar — 4.6875 x 10 -2 
Number of grid points N r — 220 
Time step At = 1.875 x 10 _1 
Propagation steps N t = 8192 


shown in panel A is that the diabatic coupling induced 
population transfer between the PESs for times t > 200. 

Our QOC goal was to maximize the nuclear probability 
flux J 2 while minimizing the flux J 3 beyond the dissocia¬ 
tive limit Rd = 8.4. To reach this aim we employed the 
following fitness function: 


J = J2 + 


1 

Js' 


(16) 


TABLE II. Boundaries of the LCP parameters space used in 
our GA approach. 


Parameter 

Minimum 

Maximum 

Eo 

0.01 

0.2 

To 

1.2 x 10 2 

9.0 x 10 2 

T 

20.0 

60.0 

UJ 

0.14 

0.16 

c 

-5.0 x 10- 7 

5.0 x 10- 7 


We analyzed the dependence of the number of LCPs, 
used in the GA optimization, on the effectiveness of the 
pulse in driving the system towards the selected target. 
To this aim, we tested three pulses built with 30, 40, and 
60 LCPs. These pulses are shown in time and frequency 
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FIG. 1. Time dependence of the populations (panel A) and 
currents (panel B) on the dissociative channels E 2 (R) and 
E 3 (R ) for the free evolution of the wavepacket. Results for 
E\(R) are not reported, as this curve is neither coupled to 
E 2 (R) nor E 3 (R). 



domains in Figures [2] and [3j Figure [3] reveals that the 
spectra of these pulses is band-width limited; this is a 
desirable feature for tailoring pulses in the laboratory. 
The spectrum shown in panel C of Figure [3] displays a 
dominant band around ujq = 0.154, frequency that is in 
the range of the transitions 3—^1 and 1 2. 


FIG. 2. Optimal pulses calculated to maximize J (Eq. 16). 
Panels A, B and C plot pulses assembled with 30, 40, and 60 
LCPs, respectively. 



FIG. 3. Power spectra of the optimal pulses plotted in Figure 

m 



Panels D, E, and F in Figure [4] present the evolution of 
the wavepacket driven by the pulse on the PES Es(R), 
E 2 {R) , and Ei(R), respectively. As shown in Figure [2b, 
the pulse displays an intense peak at t & 200, which 
drives a significant portion of the initial wavepacket to¬ 
wards the electronic ground state. This process is shown 
in more detail in panels D and F of Figure [4] Simultane¬ 
ously, the pulse induced transitions from Ei(R) to E 2 (R) 
and Es(R). 

We now analyze the efficiency of each of the pulses of Fig¬ 
ure [2] in attaining the proposed target. We perform this 
analysis in terms of time-integrated probability currents 
Ji of Eq. ( [T5| ). Panels A, B, and C of Figure [ 5 ] display the 
integrated currents at the dissociation limit (R^ = 8.4) 
on channels Ei(R), E 2 (R ), and Es(R), respectively. As 
observed, each pulse steers the system successfully to¬ 
wards the proposed goal: maximizing J 2 while minimiz¬ 
ing J3. We employed the ratio J2/J3 at time 1400 to 
quantify the efficiency of each pulse to reach the pro¬ 
posed goal. The resulting ratios J2/J3 = 10.9, 11.4, 12.3 
for pulses A, B, and C, respectively, suggest that the ef¬ 
ficiency of the optimal pulse increases with the number 
of LCPs. 

The evolution of the populations induced by the optimal 
pulse of Figure [2p is reported in Figure [6] As observed, 
at t ~ 400 the pulse transferred approximately half of 
the wavepacket from E^{R) to E\(R) and E 2 (R). In 
contrast, Figure [l]4 shows that the free evolution of the 
populations at t = 400 was barely affected by the diabatic 
coupling. 
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FIG. 4. The free evolution of the wavepacket on curves Es(R), 
E 2 (R), and Ei(R) is shown on panels A, B, and C, respec¬ 
tively. Panels D, E, and F show the time dependence of the 
wavepacket driven by the optimal pulse plotted on Figure [2p 
on curves Es(R), E 2 (R), and Ei(R), respectively. 



A. Principal component analysis 


With the aid of statistical tools we can enrich our un¬ 
derstanding of the intricate sequence of events resulting 
from the effects of the pulse and the diabatic coupling on 
the wavepacket. 

To this aim, we first built a Markov chain by record¬ 
ing the time-integrated transition moment integrals (TI- 
TMI) of each surviving individual along the GA evo¬ 
lution. Each TI-TMI is constructed as the vector 


FIG. 5. Time-integrated currents as function of the number 
of LCPs. Yellow, blue, and red lines represent the currents 
obtained with optimal pulses built with 30, 40 and 60 LCPs, 
respectively. 



FIG. 6. Time evolution of the populations on the PESs driven 
by the optimal pulse of Figure |2p. 



(Pi^2,P2^3,P3^i), where: 

Pi^j= [ T \(Xi(t)\R\Xj{t))\ 2 dt/T. ( 17 ) 

Jo 

The resulting distribution of the TI-TMIs narrows 
around the most frequent optical transitions along the 
optimal path. 

A principal component analysis (PCA) performed on this 
Markov chain provides singular values and vectors that 
can be used to enhance our understanding of the sys¬ 
tem. For instance, it supplies information about the main 
transitions induced by the pulse. Table |III| summarizes 
the PCA results for the pulse of Figure [2p when acting 
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on the molecule. The absolute value of the components 
of each singular vector measures the contribution of the 
transitions occurring along the optimal path. 


TABLE III. PCA results for the photo-processes induced 
by the pulse of Figure |2p. Pi^-j is the transition amplitude 
between PESs i and j and W 312 is the singular value for the 
corresponding process. 



Process 1 

Process 2 

Pi ->-2 

-0.150 

0.825 

E 2 —>3 

0.770 

-0.253 

ft-)-1 

-0.620 

-0.500 

W 312 

1.126 

1.100 


Table [III] lists the singular values, LF312, that quantify the 
contribution of the process along the optimal path. The 
time-integrated probability amplitudes, Pi~>j, for each 
optical transition are also reported in this table. 

We obtained two dominant processes: process 1 with sin¬ 
gular value 1.126, and process 2 with singular value 1.100. 
As observed in Table III, process 1 comprises mainly op¬ 
tical transitions P ^2 and P^i- This process is most 
likely associated with the earlier stages of the dynam¬ 
ics (£ < 500) where a significant portion of the initial 
wavepacket was driven towards the ground PES. This 
could be seen in Figure [2p. Process 2 comprises mainly 
optical transitions P^i and P\^ 2 - This process is most 
likely related with the dynamics for t > 500, where the 
pulse drives the portion of the wavepacket in the ground 
state towards its dissociation through PESs E 2 and E$. 
This could be seen in Figure [2p. 


IV. CONCLUSIONS 

The methodology proposed in this work achieves QOC 
employing GA and constraining the optimal pulse to 
analytical functions. The resulting pulse can be repro¬ 
duced in the laboratory with the current technology. 
The successful application of this QOC approach relies 
on the choice of a basis set of pulse functions. These 
functions must have spectral properties compatible with 
the transition energies of the system. 

The proposed methodology is robust and general, 
thereby allowing the use of linear superpositions of other 
experimentally available pulse functions. The suitable 
selection of the basis set of pulse functions will improve 
the performance and accuracy of the method to achieve 
more complex targets. 

The PCA performed on the Markov chain of the TI-TMIs 
shed light on the photo-dynamical processes induced 
by the optimal pulse. It proved to be an excellent tool 
for obtaining quantitative information about the optical 
processes taking place along the optimal path. 

We will apply our approach to control more complex 
systems. For instance, we will optimize sequences 
of pulses in Coherent Anti-Stokes Raman Scattering 
(CARS) j 28 * 29 * to control the dynamics in processes such 


as proton-transfer and photoisomerization. 
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